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Abstract 



Influenza viruses undergo continual antigenic evolution allowing mutant viruses to 

evade immunity acquired by the host population to previous virus strains. Antigenic 

phenotype is often assessed through pairwise measurement of cross-reactivity between 

(y-) influenza strains using the hemagglutination inhibition (HI) assay. Here, we extend 



previous approaches to antigenic cartography, which seeks to place strains on an anti- 
genic map, such that distances on this map best recapitulate titers observed across 
(f~) multiple HI assays. In our model, we simultaneously characterize antigenic and ge- 

netic evolution by including an evolutionary model in which antigenic location diffuses 
over a shared virus phylogeny. Using HI data for four lineages of influenza, encom- 
•*~j passing influenza A subtypes H3N2 and H1N1, and influenza B lineages Victoria and 

rN Yamagata, we determine average rates of antigenic drift for each lineage, as well as 

year-to-year variability in the rate of drift. Through comparison with epidemiologi- 
cal data, we demonstrate a year-to-year correlation between drift and incidence and 
present evidence that antigenic drift mediates interference between influenza lineages. 
We investigate the selective underpinnings for differing antigenic dynamics across lin- 
eages and show that A/H3N2 benefits from both a higher influx of new antigenic 
mutations and also from more efficient conversion of antigenic variation into fixed 
differences. This work does much to elucidate the antigenic dynamics of influenza 
lineages, but also allows for substantial future advances in investigating the dynam- 
ics of influenza and other antigenically- variable pathogens by providing a model that 
intimately combines molecular and antigenic evolution. 



Introduction 

Seasonal influenza infects between 10% and 20% of the human population every year, 
causing an estimated 250,000 to 500,000 deaths annually [lj. While individuals develop 
long-lasting immunity to particular influenza strains after infection, antigenic mutations to 
the influenza virus genome result in proteins that are recognized to a lesser degree by the 
human immune system, leaving individuals susceptible to future infection. The influenza 
virus population continually evolves in antigenic phenotype in a process known as antigenic 
drift. A large proportion of the disease burden of influenza stems from antigenic drift; it is 
why vaccines remain only transiently effective. A thorough understanding of the process 
of antigenic drift is essential to public health efforts to control mortality and morbidity 
through the use of a seasonal influenza vaccine. 

Before 2009, there were four major clades or lineages of influenza circulating within the 
human population: the H3N2 and H1N1 subtypes of influenza A, and the Victoria and Ya- 
magata lineages of influenza B. In the case of influenza A, subtypes A/H3N2 and A/H1N1 
refer to the genes, hemagglutinin (H or HA) and neuraminidase (N or NA), that are 
primarily responsible for the antigenic character of a strain. In the case of influenza B, 
Victoria (B/Vic) and Yamagata (B/Yam) refer to antigenically distinct lineages which 
diverged from a single lineage prior to 1980 |2|. Mutations to the HAl region of the 
hemagglutinin protein are thought to drive the majority of antigenic drift in the influenza 
virus [3l[4J . Experimental characterization of antigenic phenotype is possible through the 
hemagglutination inhibition (HI) assay |5|, which measures the cross- reactivity of one 
virus strain to serum raised against another strain through challenge or vaccination. Sera 
from older strains react poorly with more evolved viruses resulting in new strains having 
a transmission advantage over previously established strains. 

The results of many HI assays across a multitude of virus strains of a single subtype 
can be combined to yield a two-dimensional map, quantifying antigenic similarity and 
distance p] . The antigenic map of influenza A/H3N2 has shown largely linear movement 
of the influenza virus population since its introduction in 1968. Evolution of antigenic 
phenotype appears punctuated with episodes of more rapid innovation interspersed by 
periods of stasis, while genetic evolution appears more continuous |6|, suggesting that 
a relatively small number of genetic changes or combinations of genetic changes may 
drive changes in antigenic phenotype. The process of antigenic drift results in the rapid 
turnover of the virus population. Although mutation occurs rapidly, genetic diversity 
among contemporaneous viruses is low and phylogenetic analysis shows a characteristically 
'spindly' tree with a single predominant trunk lineage and transitory side branches that 
persist for only 1-5 years [71. 

Previously, the antigenic and genetic patterns of influenza evolution have been analyzed 
essentially in isolation. An antigenic map is constructed from a panel of HI measurements, 
and a phylogenetic tree is constructed from sequence data. However, the opportunity for 
a combined approach exists as both the antigenic map and the phylogenetic tree often 
contain many of the same isolates. Here, we implement a flexible Bayesian approach to 
characterize jointly the antigenic and genetic evolution of the influenza virus population. 
We apply this approach to investigate the dynamics of A/H3N2, A/H1N1, B/Vic and 



B/Yam viruses, and, for the first time, present detailed reconstructions of the antigenic 
dynamics of all four circulating influenza clades. 



Results and discussion 

Antigenic and evolutionary cartography 

In order to assess patterns of antigenic evolution among influenza strains, we imple- 
mented a Bayesian probabilistic analog of multidimensional scaling (MDS), referred to 
here as BMDS (see Materials and methods). In this model, viruses and sera are given 
iV-dimensional locations, thus specifying an 'antigenic map', such that distances between 
viruses and sera in this space are inversely proportional to cross-reactivity. In the BMDS 
model, a map distance of one antigenic unit translates to an expectation of a 2-fold drop 
in HI titer between virus and sera. Maps that produce pairwise distances most congruent 
with the observed titers will have a high likelihood and will be favored by the BMDS 
model. We integrate over sources of uncertainty, such as antigenic locations, in a flexible 
Bayesian fashion. We apply this model to HI measurements of virus isolates against post- 
infection ferret antisera for influenza A/H3N2, A/H1N1, B/Vic and B/Yam (see Materials 
and methods). 

We begin with Bayesian analogs of the models used by Smith et al. J6], in which viruses 
and sera are represented as iV-dimensional locations, expected titer is relative to the 
maximum titer of a particular ferret serum and virus and serum locations follow a diffuse 
prior. After comparing models of differing dimensions, Smith et al. p] arrive at a 2D 
model as the preferred model for their data. Smith et al. p\ implement a form of MDS, 
seeking to optimize virus and serum locations such that the sum of squared errors between 
expected and observed titers is minimized (eq. k^ of Materials and methods). Here, in 
implementing BMDS, we provide a likelihood function for the probability of observing HI 
data given virus and serum locations (eq. [8] of Materials and methods) and seek to estimate 
model parameters through Bayesian inference using Markov chain Monte Carlo (MCMC). 
However, the basic antigenic model describing drop in HI titer as proportional to Euclidean 
distance between virus and serum locations is identical between these methods. 

We test model performance by constructing training datasets representing 90% of the HI 
measurements for each of the four influenza lineages and test datasets representing the 
remaining 10% of the measurements for each lineage. By fitting the BMDS model to 
the training dataset, we are able to predict HI titers in the test dataset and compare 
these predicted titers to observed titers. We find that a two dimensional model has better 
predictive power than models of lower of higher dimension in all four influenza lineages 
(models 1-5; Table IT]). We find that this 2D model performs well, yielding an average 
absolute predictive error of between 0.78 and 0.91 log2 HI titers across influenza lineages 
(model 2; Table [T]), in line with the results of Smith et al. (6J. Consequently, we specify 
a two dimensional model in all subsequent analyses. The finding of a low-dimensional 
map across influenza lineages extends previous studies in A/H3N2 |6| and remains an 
interesting and fundamental empirical observation. 



Table 1. Average absolute prediction error of log2 HI titer for test data across models 
and datasets. 















Test error 




Model 


Dimen 


Location prior 


Serum effects 


Virus effects 


A/H3N2 


A/H1N1 


B/Vic 


B/Yam 


1 


ID 


Uninformed 


Fixed 


None 


1.35 


0.94 


0.90 


1.08 


2 


2D 


Uninformed 


Fixed 


None 


0.91 


0.78 


0.82 


0.90 


3 


3D 


Uninformed 


Fixed 


None 


0.93 


0.80 


0.85 


0.92 


4 


4D 


Uninformed 


Fixed 


None 


0.98 


0.84 


0.90 


0.97 


5 


5D 


Uninformed 


Fixed 


None 


1.04 


0.89 


0.98 


1.04 


6 


2D 


Drift 


Fixed 


None 


0.91 


0.75 


0.77 


0.83 


7 


2D 


Diffusion/Drift 


Fixed 


None 


0.89 


0.74 


0.74 


0.83 


8 


2D 


Diffusion/Drift 


Estimated 


None 


0.77 


0.73 


0.66 


0.75 


9 


2D 


Diffusion/Drift 


Estimated 


Estimated 


0.76 


0.71 


0.64 


0.72 



Previous work on influenza antigenic and genetic evolution has examined antigenic and 
genetic patterns in isolation, and then compared the results to assess congruence of these 
two aspects of evolution [6,8,9|. Here, we simultaneously model antigenic and genetic 
evolution by adopting an evolutionary diffusion process [10], wherein a virus's antigenic 
character state evolves along branches of the phylogenetic tree according to a Brownian 
motion process (see Materials and methods). The phylogenetic diffusion process acts 
as a prior on virus locations, so that genetically similar viruses are expected to share 
similar antigenic locations. The antigenic diffusion process includes both systematic drift 
with time and covariance induced by phylogenetic proximity. We examine the effects of 
including only systematic drift (model 6; Table [l]) and systematic drift plus phylogenetic 
diffusion (model 7; TablefTl), finding a small increase in predictive accuracy of between 0.02 
and 0.08 log2 HI titers when both processes are included. The systematic drift process 
informs virus and serum locations by dates of isolation and the phylogenetic diffusion 
process informs virus locations by genetic sequences. Thus, in these models, antigenic 
locations are inferred using both genetic data and HI data and will differ from locations 
inferred from HI data alone. If HI data is rich, then we expect minor differences in antigenic 
locations with the inclusion of genetic data (as may be the case for A/H3N2), while if HI 
data is spare, then we expect genetic data to play a larger role in determining antigenic 
locations (as may be the case for B/Vic and B/Yam). 

We further extend the model by estimating the strength of overall reactivity of each serum 
rather than fixing this at its maximum titer, and additionally, by estimating the strength 
of reactivity of each virus isolate. We refer to these estimates as serum effects and virus 
effects, respectively. We find that including these effects decreases test error further, 
ranging in improvement from 0.03 to 0.13 log2 HI titers (models 8 and 9; TablefTl). 

We find that the average absolute error in predicted log2 HI titer is nearly constant with 
antigenic distance (Figure fljA.), thus supporting our model assumption that the drop in log2 
titer is proportional to the Euclidean distance separating viruses and sera on the antigenic 
map. Additionally, we find that the absolute error in predicted titer is nearly constant with 
time (Figure [Ip). Antigenic locations inferred by the model are well resolved; estimates 
of antigenic distance between pairs of viruses show relatively little variation across the 



posterior. We estimate that virus distances have, on average, a 50% credible interval of 
±0.45 antigenic units for A/H3N2, ±0.57 units for A/H1N1, ±0.76 units for B/Vic, and 
±0.65 units for B/Yam. 
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Figure 1. Absolute error in log2 HI titer against antigenic distance and date of virus 
sampling in influenza A/H3N2. (A) Points show the absolute deviation from the predicted 
log2 HI titer of the antigenic model, computed from virus location Xj, serum location y^, virus 
effect m and serum effect Sj, based upon antigenic distance between virus and serum Sij. The red 
line shows a LOESS regression curve of these points. (B) Points show the absolute deviation from 
the predicted log2 HI titers of the antigenic model, based upon the date of virus sampling. The 
red line shows average deviations for each year. Relationships for A/H1N1, B/Vic and B/Yam are 
similar. 



We find strong correspondence between our results and previous results by Smith et al. p] , 
with equivalent models producing globally consistent antigenic maps and other models 
producing locally consistent maps with a small degree of global inconsistency (Figures pl- 
11 ). Differences in the MDS and BMDS approaches reflect greater philosophical differences 
between maximum-likelihood and Bayesian statistical approaches, with the former seeking 
the single most likely explanation for the data, and the latter seeking to fully characterize 
model uncertainty. If a single map is desired, then the MDS method may be preferred, 
while the BMDS method provides flexibility, allowing other sources of information, such 
as genetic data, to be more easily incorporated. 



Antigenic evolution across influenza lineages 



Through our analysis, we reveal the antigenic, as well as evolutionary, relationships among 
viruses in influenza A/H3N2, A/H1N1, B/Vic and B/Yam, quantifying both antigenic 
and evolutionary distances between strains (Figure^]). Over the time period of 1968 to 
2011, influenza A/H3N2 shows substantially more antigenic evolution than is exhibited 
by A/H1N1 over the course of 1977 to 2009 or B/Vic and B/Yam over the course of 
1986 to 2011. We observe prominent antigenic clusters in A/H3N2 and A/H1N1, but less 
prominent, though still apparent, clustering in B/Vic and B/Yam. Antigenic clusters show 
high genetic similarity, so that we observe very few mutation events leading to each cluster, 
rather than the repeated emergence of clusters. This analysis makes the fate of antigenic 
clusters obvious, with two clusters in A/H3N2 (Victoria/75 and Beijing/89) appearing to 
be evolutionary dead-ends. Labeling of prominent antigenic clusters in figures [2] and [3] 



is intended as a rough guide for orientation and not as exhaustive catalog of antigenic 
variation. 
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Figure 2. Antigenic locations of A/H3N2, A/H1N1, B/Vic and B/Yam viruses show- 
ing evolutionary relationships between virus samples. Circles represent a posterior sam- 
ple of virus locations and have been shaded based on year of isolation. Antigenic units rep- 
resent two-fold dilutions of the HI assay. Absolute positioning of lineages, e.g. A/H3N2 and 
A/H1N1, is arbitrary. Lines represent mean posterior diffusion paths when virus locations are 
hxed. Prominent antigenic clusters are labeled after vaccine strains present within clusters, and 
are abbreviated from Hong Kong/68, England/72, Victoria/75, Bangkok/79, Sichuan/87, Bei- 
jing/89, Beijing/92, Wuhan/95, Sydney/97, Fujian/02, California/04, Wisconsin/05, Brisbane/07, 
Perth/09 (A/H3N2), USSR/77, Singapore/86, Beijing/95, New Caledonia/99, Solomon Islands/06 
(H1N1), Victoria/87, Hong Kong/01, Malaysia/04, Brisbane/08 (Vic), Yamagata/88, Shanghai/02, 
Florida/06, Wisconsin/10 (Yam). 



HI assays lack sensitivity beyond a certain point, so that for A/H3N2, cross-reactive mea- 
surements only exist between strains sampled at most 14 years apart, leaving only thresh- 
old titers, e.g. '<40', in more temporally distant comparisons. Because of the threshold 
of sensitivity of the HI assay, it's difficult to distinguish a linear trajectory in 2D anti- 
genic space from a slightly curved trajectory. To solve this problem of identifiability, we 
assumed a weak prior that favors linear movement in the 2D antigenic space (present in 
models 6 through 9; Table [l]), with the slope of the linear relationship and the precision 
of the relationship incorporated into the Bayesian model (see Materials and methods). 
Because of this, we interpret map locations locally rather than globally, and assess rates 
of antigenic movement without making strong statements about the larger configuration 
under which the movement occurs. 

We find that influenza A/H3N2 evolved along antigenic dimension 1 at a rate of 1.01 
antigenic units per year, with a 95% highest posterior density (HPD) interval of 0.98-1.04 
(Figure [3]). However, we observe occasional large jumps in antigenic phenotype (Figure [3]), 
corresponding to cluster transitions identified by Smith et al. j6j. Most variation is con- 
tained within the first antigenic dimension, but dimension 2 occasionally shows variation 
when two antigenically distinct lineages emerge and transiently coexist (Figure [2]) , as is 
the case with the previously identified Beijing/89 and Beijing/92 clusters. 
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Figure 3. Antigenic drift of A/H3N2, A/H1N1, B/Vic and B/Yam viruses showing 
evolutionary relationships between virus samples. Antigenic drift is shown in terms of 
change of location in the first antigenic dimension through time. Circles represent a posterior sam- 
ple of virus locations and have been shaded based on year of isolation. Antigenic units represent 
two- fold dilutions of the HI assay. Relative positioning of lineages, e.g. A/H3N2 and A/H1N1, in 
the vertical axis is arbitrary. Lines represent mean posterior diffusion paths when virus locations 
are fixed. Prominent antigenic clusters are labeled after vaccine strains present within clusters, 
and are abbreviated from Hong Kong/68, England/72, Victoria/75, Bangkok/79, Sichuan/87, Bei- 
jing/89, Beijing/92, Wuhan/95, Sydney/97, Fujian/02, California/04, Wisconsin/05, Brisbane/07, 
Perth/09 (A/H3N2), USSR/77, Singapore/86, Beijing/95, New Caledonia/99, Solomon Islands/06 
(H1N1), Victoria/87, Hong Kong/01, Malaysia/04, Brisbane/08 (Vic), Yamagata/88, Shanghai/02, 
Florida/06, Wisconsin/10 (Yam). 



We find that other clades of influenza evolved in antigenic phenotype substantially slower 
than A/H3N2 (Figure pjl). Influenza A/H1N1 evolved at a rate of 0.62 units per year 
(HPD 0.56-0.67), but showing a similar pattern of punctuated antigenic evolution with 
occasional larger jumps in phenotype, such as the emergence of the Solomon Islands/06 
cluster. Influenza B/Victoria also evolved relatively slowly, with an average rates of 0.42 
(HPD 0.32-0.51) units per year. Influenza B/Yamagata evolved slower still, with an 



average rate of 0.32 (HPD 0.25-0.39) units per year. Punctuated evolution is less obvious 
in B/Yam and B/Vic compared to A/H3N2 and A/H1N1, but antigenic clusters are still 
apparent, with recent transitions to the Brisbane/08 cluster in B/Vic fiT] and to the 
Wisconsin/10 cluster in B/Yam [12] . Interestingly, a minor lineage of B/Vic, denoted 
B/Hubei-Songzi/51/2008 11], has persisted through 2011, while remaining antigenically 
distinct from B/Brisbane/60/2008 viruses (Figure pi). 

We sought to summarize year-to-year patterns of antigenic drift by calculating the dif- 
ference in mean location between consecutive years (Figure H| . We estimate year-to-year 
antigenic drift for years 1992 to 2011 by calculating the average location along dimension 
1 of phylogenetic lineages present in the tree at year i and comparing this location to the 
average location of phylogenetic lineages present in the tree at year i — 1. There may often 
be large discontinuities in virus locations across the population; our use of difference in 
mean location is meant to capture both the distance between antigenic clusters and also 
the change in cluster frequency over consecutive years. We observe greater heterogeneity 
in year-to-year antigenic drift in type A than in type B lineages (Figure El) , with stan- 
dard deviation of year-to-year antigenic drift equal to 0.97 units in A/H3N2, 0.66 units in 
A/H1N1, 0.46 units in B/Vic and 0.26 units in B/Yam. This analysis classifies drift only 
to the level of consecutive years; some coarse-graining of the timings of transition events 
will necessarily occur. 



Epidemiological consequences of antigenic drift 

We investigate the relationship between rates of antigenic drift and USA incidence in 
A/H3N2, A/H1N1, B/Vic and B/Yam. We take estimates of the rate of antigenic drift 
from the preceding section and compare these rates to measurements of relative incidence 
for each influenza clade calculated from the proportion of influenza viruses attributable 
to each clade (see Materials and methods). We analyze incidence from the 1998/1999 to 
the 2008/2009 seasons to avoid possible complications from the 2009 pandemic, finding 
that A/H3N2 accounted for 56%, A/H1N1 for 20%, B/Vic for 10% and B/Yam for 14% 
of viruses detected in the USA during this time period. We find a significant correla- 
tion between rate of antigenic drift and relative incidence across the four clades (Pearson 
correlation, r = 0.96, p = 0.021). 

We follow-up this analysis with a more detailed analysis of year-to-year variation in anti- 
genic drift and clade-specific incidence in the USA. We measure year-to-year antigenic drift 
within each clade as described in the previous section and shown in Figure [4] We measure 
seasonal incidence of each clade by taking average influenza-like illness (ILI) percentage 
and multiplying this by the proportion of viruses attributable to a clade for each season. 
This measure of incidence has previously been shown to have have predictive power in the 
analysis of seasonal influenza trends [13]. In correlating antigenic drift to seasonal inci- 
dence we standardize this measure of incidence across seasons and across clades to measure 
effects in terms of standard deviations of incidence. We find that years with high levels of 
antigenic drift coming into an influenza season, e.g. drift in A/H3N2 from 2000 to 2001, 
show high levels of incidence in that season, e.g. incidence of A/H3N2 in the 2001/2002 
season (FigureEl FigureplA.). Here, within-clade drift explains 35% of the variance in stan- 
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Figure 4. Year-to-year antigenic drift between 1992 and 2009 in A/H3N2, A/H1N1, 
B/Vic and B/Yam viruses. (A) Time series of year-to-year antigenic drift across influenza 
clades from 1992 to 2011 and seasonal incidence from 1998/1999 to 2008/2009. Colored lines 
represent year-to-year antigenic drift, and drift for a virus clade for year i is measured as the 
mean of antigenic dimension 1 of phylogenetic lineages in year i compared to the mean of antigenic 
dimension 1 of phylogenetic lineages from the previous year % — 1. For example, 2000 represents 
difference in antigenic dimension 1 between 1999 and 2000. Error bars represent 50% Bayesian 
credible intervals of year-to-year drift. Gray dotted lines represent clade-specihe USA incidence 
taken as average influenza-like illness (ILI) multiplied by proportion of viruses attributable to a 
clade for each season. Here, 2000 represents the 2000/2001 influenza season. (B) Histograms of 
year-to-year antigenic drift across influenza clades summarizing observations from 1992 to 2011. 



dardized incidence and represents a significant correlation (Pearson correlation, r 
p = 1.6x HT 5 ). 
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We further show that incidence within a clade decreases with antigenic drift in other clades, 
e.g. drift in A/H1N1, B/Vic and B/Yam from 2000 to 2001 decrease incidence of A/H3N2 
in the 2001/2002 season (Figure [5t$)- This also represents a significant relationship (Pear- 
son correlation, r = —0.20, p = 0.024). We fit a linear model that predicts standardized 
incidence based upon drift within a clade, drift in other clades of the same virus type, and 
drift in clades of differing type (Table ^1). For example, standardized incidence in A/H3N2 
would be explained by drift within A/H3N2, drift in A/H1N1 and by drift across types 
in B/Vic and B/Yam. We compared adjusted R 2 across models with different regression 
coefficients included, and arrived at a best-fitting model that explains 46% of the variance 
in standardized incidence and includes the terms (3 C , /3 S and fit, but does not include an 
intercept term. We find that 1 unit of antigenic drift within a clade increases incidence by 
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Figure 5. Relationship between antigenic drift within and between clades and seasonal 
USA incidence for years 1998 to 2009. (A) Standardized incidence within a clade compared 
to antigenic drift within the same clade. For example, incidence of A/H3N2 in season 1998/1999 
is measured against year-to-year antigenic drift of A/H3N2 from 1997 to 1998. (B) Standardized 
incidence within a clade compared to antigenic drift in other clades. For example, incidence of 
A/H3N2 in season 1998/1999 is measured against year-to-year antigenic drift of A/H1N1 from 
1997 to 1998, as well as antigenic drift of B/Vic from 1997 to 1998. Standardized incidence is 
measured as yearly ILI percentage multiplied by proportion of viruses identified in a particular 
clade, and then standardized so that mean and variance across seasons and clades is equal to and 
1 respectively. Colors represent the clade in which incidence is measured; light blue for A/H3N2, 
purple for A/H1N1, orange for B/Vic and red for B/Yam. 

1.26 standard deviations, while 1 unit of drift in a clade of the same type decreases inci- 
dence by 0.61 standard deviations and 1 unit of drift in clades of differing type decreases 
incidence by 0.32 standard deviations. The term for across lineage interference fit is small 
in magnitude, but remains statistically significant (Table ^1). 

Table 2. Estimates of regression coefficients for the relationship of clade-specific USA 
incidence to antigenic drift within and across clades, where standardized incidence in 
A/H3N2 is predicted by fi c xh3 + fi s #H1 + fit ^Vic + fit #Yam an d standardized incidence in 
other clades follows an analogous relationship. 



Parameter Estimate Standard error p-value 



fic 


1.26 


0.22 


1.1 x 10~ 6 


fis 


-0.61 


0.22 


0.008 


fit 


-0.32 


0.11 


0.005 



These findings clearly demonstrate the epidemiological importance of antigenic drift and 
show that incidence across A/H3N2, A/H1N1, B/Vic and B/Yam is strongly influenced 
by levels of antigenic drift within each clade, consistent with antigenic drift increasing 
the susceptible pool available to a clade. The finding of a negative correlation between 
incidence within one clade and antigenic drift in other clades supports the hypothesis 
that epidemiological interference exists between clades [13[|14| . Interestingly, we find that 
interference between sister clades, e.g. B/Vic and B/Yam, appears stronger than interfer- 
ence between more divergent clades, e.g. B/Vic and A/H3N2. If interference is mediated 
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by CTL epitopes, then we would expect interference to decrease with evolutionary dis- 
tance (Tsl. 



Furthermore, these findings suggest that the predominant clade of a particular influenza 
season could be estimated ahead of time by examining the degree of variation in antigenic 
drift between influenza clades. However, this prediction would not be perfect, and 56% 
of the variance in incidence remains unexplained by these estimates of antigenic drift. 
Some of this variance may be explained by incorporating geographic variation in year-to- 
year antigenic drift; although there exists significant geographic variation in clade-specific 
incidence [16] , we calculated antigenic drift from global samples, rather than USA-specific 
samples. 

Conversion of antigenic mutation into fixed differences 

We observe a faster rate of antigenic drift in influenza A/H3N2 than in A/H1N1 or in either 
lineage of influenza B. Previous work using general epidemiological models has suggested 
that rates of antigenic drift may be influenced by both the fundamental reproductive num- 



ber Rq and the rate at which mutation decreases cross- immunity 117 , 18 . Owing to these 
and similar findings, subsequent work ascribed differences in the evolution of A/H3N2 
relative to A/H1N1 and influenza B to either greater Ro or greater mutation rate [14||19] . 
Here, we examine the conversion of recently acquired antigenic mutation into fixed differ- 
ences by comparing rates of antigenic diffusion to the degree of persistence on phylogeny 
branches in A/H3N2, A/H1N1, B/Vic and B/Yam (Table pi). To facilitate comparisons, 
we measure the diffusion coefficient D |20| characterizing the antigenic movement of each 
influenza lineage (see Materials and methods). This coefficient is measured in terms of 
units 2 /year and can be used to determine the diffusion length scale, that is the expected 
relationship between time and distance traveled, which is 2\ftD for a two-dimensional 
diffusion. Higher diffusion coefficients mean a lineage is diffusing faster through antigenic 
space. We compare diffusion coefficients between trunk branches and side branches, where 
trunk branches are defined as branches ancestral to the most recent virus sample and side 
branches are defined as those branches not ancestral to the most recent virus sample. 

Table 3. Estimates of fundamental diffusion coefficient D across influenza lineages and 
phylogenetic branches, including posterior means and 50% highest posterior density 
intervals. 

Lineage Trunk branch D Side branch D Ratio 

A/H3N2 0.90 (0.79-0.97) 0.62 (0.57-0.65) 1.45 (1.31-1.57) 

A/H1N1 0.62 (0.50-0.71) 0.44 (0.38-0.48) 1.45 (1.18-1.70) 

B/Vic 0.63 (0.51-0.71) 0.62 (0.54-0.68) 1.02 (0.87-1.14) 

B/Yam 0.35 (0.29-0.40) 0.33 (0.28-0.36) 1.09 (0.93-1.23) 

Mutations must proceed through the sieve of natural selection, first appearing as transient 
differences restricted to side branches, while advantageous mutations will preferentially 
fix in the virus population and thus accrue on the trunk of the phylogeny. We expect 
the diffusion coefficient on side branches to more closely resemble the underlying rate of 
antigenic mutation, before selection can strongly influence observed rates. Here, we find 
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that A/H3N2 and B/Vic have similar diffusion coefficients on side branches, which are 
larger than the diffusion coefficients of A/H1N1 and B/Yam (Table pi). Thus, we expect 
that new antigenic mutations may occur at a faster rate in A/H3N2 and B/Vic than in 
A/H1N1 or B/Yam. We find that trunk branches in both A/H3N2 and A/H1N1 evolve in 
antigenic phenotype faster than side branches (Table [3]) , suggesting that positive selection 
is promoting antigenic movement in these lineages. Interestingly, B/Vic and B/Yam show 
less of a signature of positive selection, with similar diffusion coefficients on trunk and 
side branches (Table pjj). Thus, although drift forward along dimension 1 is driven by 
selection in all lineages (Figure [3]), raw antigenic diffusion appears to be selected for more 
strongly in influenza A than in influenza B, giving influenza A an advantage in overall 
rate of antigenic drift. Additionally, it appears that A/H3N2 and B/Vic benefit from a 
greater influx of new antigenic mutations, giving them a similar advantage over their sister 
lineages A/H1N1 and B/Yam, respectively. 

In this context, it is especially interesting that adaptive evolution in the HA1 protein deter- 
mined through analysis of nonsynonymous and synonymous substitutions is substantially 
stronger in A/H3N2 than in A/H1N1 [21] ; H3 appears to accumulate adaptive fixations 
at ~150% the rate of HI [22]. Here, we observe that A/H3N2 has a diffusion coefficient 
~141% greater than A/H1N1 along side branches of the phylogeny. These results suggest 
that the greater rate of adaptive fixation in A/H3N2 may be simply due to a greater 
abundance of new antigenic mutations to fix. 

We suggest that the larger signal of positive selection for antigenic diffusion in influenza 
A relative to influenza B could be due to either a greater fundamental reproductive num- 
ber Rq or to pleiotropic effects associated with antigenic change. In the latter case, we 
expect that many antigenic variants may be at an inherent disadvantage against other 



strains from decreased protein function |23 24 or other effects, even though they benefit 



from a transmission advantage mediated by host immunity. These mutations may not 
be promoted by natural selection to the same extent as antigenic variants that maintain 
functional equivalence. Additionally, other work has shown strong pleiotropic constraints 
in the evolution of oseltamivir resistance in A/H1N1 [25] . Still, it is clear that more work 
is necessary to fully clarify the mechanistic underpinnings of the relative evolutionary 
success of different influenza lineages. 



Conclusions 

In this study we provide a foundation for evolutionary antigenic cartography, which seeks 
to simultaneously assess antigenic phenotype and antigenic evolution. We use this ap- 
proach to characterize competitive dynamics across influenza clades A/H3N2, A/H1N1, 
B/Vic and B/Yam and show that antigenic evolution drives strain replacement in each 
clade and interference between clades. We find that varying levels of mutational input and 
varying efficiencies of the conversion of antigenic polymorphism into fixed differences re- 
sults in different rates of antigenic drift across influenza clades. Influenza A/H3N2 benefits 
from a higher rate of new antigenic mutation coupled with strong adaptive evolution and 
evolves in antigenic phenotype faster than A/H1N1, B/Vic and B/Yam. Correspondingly, 
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we observe substantially greater levels of incidence in A/H3N2 than in other influenza 
clades. We suggest that antigenic evolution strongly influences competitive dynamics 
both within and between influenza clades. 

The statistical framework presented here represents a baseline to which further advance- 
ments in modeling antigenic phenotype and evolution may be made. For example, our 
likelihood-based model facilitates the inclusion of possible covariates affecting immuno- 
logical titer, which could include experimental factors such as red blood cell type used 
in the HI assay [26] and whether oseltamivir is included in the HI reaction [27]. Addi- 
tionally, this framework should be ideally suited to uncovering genetic determinants of 
antigenic change, as both the sequence state and antigenic location of internal nodes in 
the phylogeny may be estimated. In this fashion, it should be possible to correlate se- 
quence substitutions directly to antigenic diffusion. Here, methods to efficiently calculate 
branch- and site-specific codon substitutions may be necessary to easily correlate sequence 
substitutions with antigenic change |28|. 

A deeper understanding of the process of antigenic evolution in the human influenza virus 
may eventually prove to be crucial in improving vaccine strain selection. 



Materials and methods 

Antigenic cartography 

Antigenic characteristics of viral strains are often assessed through immunological assays 
such as the hemagglutination inhibition (HI) assay (5j. At heart, these assays compare 
the reactivity of one virus strain to antibodies raised against another virus strain via 
challenge or vaccination. In the case of HI, the measurement of cross-reactivity takes the 
form of a titer representing the dilution factor at which serum raised against a particular 
virus ceases to be effective at inhibiting the binding of another virus to red blood cells. 
These factors are commonly assessed by serial dilution, so that HI titers will form a log 
series, 40, 80, 160, etc .... Because experimental HI titers typically differ by factors of 
two, we find it convenient to work in log2 space and represent the titer of virus i against 
serum j as H^ = log 2 (HI titer), i.e. a titer of 160 has H^ = 7.32. Due to experimental 
constraints, most comparisons cannot be made, leading to a sparse observation matrix 
H = {Hij}. Further, measurements are usually interval and truncated, e.g. inhibition 
may cease somewhere between the serial titers of 160 and 320, or inhibition may be absent 
at all titers assayed, suggesting a threshold somewhere between and 40. 



Previous work |6l|29] has used multidimensional scaling to place viruses and sera on an 
'antigenic map'. These methods heuristically optimize locations of viruses and sera by 
seeking to minimize the sum of squared errors between titers predicted by map locations 
and observed titers. Antigenic maps produced by these methods have proved useful in 
categorizing virus phenotypes [6] , but the extension of these methods to integrate genetic 
data remains notably lacking. 

Here, we follow previous models in representing antigenic locations as points in a low P- 
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dimensional antigenic map. One of our initial goals is to find an optimal projection of the 
high-dimensional distance matrix H into this lower dimensional space. We conduct this 
projection using Bayesian multidimensional scaling (BMDS) [ 30 1 in which we construct a 
probabilistic model to quantify the fit of a particular configuration of cartographic locations 
to the observed matrix of serological measurements. Typically, P = 2, but higher or lower 
dimensions may better reflect the data. 

Let Xj E K represent the cartographic location of virus i for i = 1, . . . , n, so that Xj = 
(xn,Xi2)' for P = 2. Similarly, let y,,- represent the cartographic location of serum j 
for j = l,...,k, so that yj = (yji,yj2)' for P = 2. For notational compactness, we 
collect together all virus coordinates into an nx P matrix X = (xi, . . . , x n )' and all serum 
coordinates into anfcxP matrix Y = (yi, . . . , y^)'. Virus and serum may be isolated from 
/ raised against the same strain and have different cartographic locations, and separate 
serum isolates raised against the same strain may also have different cartographic locations. 
This gives a set of distances between virus and serum cartographic locations 

Sij = ||xj — y j 1 1 2 , (1) 

where || • H2 is an L2 norm. 

Traditional approaches to antigenic cartography J6j begin by defining immunological dis- 
tance as 

y — ^j y ' \ ) 

where H^ is the log2 titer of virus i against serum j and serum effect Sj = m&x(Hij, . . . , H n j) 
is fixed. In following multidimensional scaling (MDS), these approaches attempt to opti- 
mize over unknown X and Y such that 



(«)6I 



i] dij) (3) 



is minimized, where I = {(i, j) : H^ is measured}. In the case of threshold measurements, 
this error function is modified slightly; see [6] for further details. 

Here, we instead assume a probabilistic interpretation in which an observed titer is nor- 
mally distributed around its cartographic expectation with variance (p 2 , 

Hij^Nisj-Sij,^ 2 ). (4) 

Consequently, the likelihood of observing an exact titer given the placement of antigenic 
locations is 

/,W = *( g * + ^~ Ji ), (5) 

where </>(•) represents the standard normal probability density function (PDF). Previous 
BMDS has employed a sampling density truncated to strictly positive quantities since 
dij are directly observed, non-negative quantities. In the antigenic setting, these remain 
random and can be negative since neither sj is known nor is Hij observed with much 
precision. 
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Figure 6. Likelihood of HI titers in the BMDS model. Here we show the likelihoods 



of observing three different outcomes given Si 



f 



0.95 and 



log 2 (1280) = 10.32. 



The likelihood of observing a threshold titer of '<40' is equal to the lower tail of the probability 
density function /_j(5.32) = 0.146. The likelihood of observing a point measurement with an exact 
inhibiting titer of '90.5' is equal to the density function /| (6.5) = 0.413. The likelihood of observing 
an interval measurement with an inhibiting titer somewhere between '160' and '320' is equal to 
/u(7.32) = 0.129. 

HI assays sometimes show no inhibition at all measured titrations, e.g. a measurement can 
be reported as '<40'. In this case, the likelihood of observing the threshold measurement 
follows the cumulative density of the lower tail of the normal distribution 



/ J (fly)=*('^ ± ^ 



(6) 



where $(•) represents the standard normal cumulative distribution function (CDF). Al- 
though it is simplest to assume that immunological measurements represent point esti- 
mates, it seems more natural to assume that the threshold for inhibition occurs between 
two titers, e.g. we observe inhibition at 1:160 dilution and no inhibition at 1:320 dilution. 
Rather than taking the HI titer as 160, we can instead treat this as an interval measure- 
ment, assuming that the exact titer for inhibition would occur somewhere between 160 
and 320. HI titers are usually reported as the highest titer that successfully inhibits virus 
binding, so that in this case, we calculate the likelihood of an interval measurement as 



fu(H., 



i.i > 



$ 



Ha + Si 



Sj + 1 
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Ha + 5. 
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(7) 



These likelihoods are illustrated in Figure [6j Throughout our analyses, we use interval 
likelihoods /u rather than point likelihoods f\ unless otherwise noted. 

We calculate the overall likelihood by multiplying probabilities of individual measurements 



L(X,Y)= \\ f(K 



»J7' 



(8) 



(*J)6I 
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using probability functions /| , /j and / u as appropriate. We begin by assuming indepen- 
dent, diffuse normal priors on virus and serum locations 

Xj ~ Normal(/i, S) 

yj ~ Normal(/x,£), (9) 

where // = (0, . . . , 0)' and 5] is a diagonal matrix with diagonal elements all equal to 10000. 

Virus and serum effects 

The preceding model represents immunological distance as a drop in titer against the most 
reactive comparison for a particular serum. However, this model may be biased in some 
circumstances. In one example, if a particular serum j is only measured against distant 
viruses, its maximum titer will be artificially low and the likelihoods concerning this serum 
will appear poor. To address this issue, we relax the assumption of fixed Sj values and treat 
the expected log2 titer when 5ij = as a random variable. In this case, Hij still follows 
equation El with expectation Sj — 5ij, but the vector of 'serum effects' s = (si,...,Sfc) 
is random and estimated rather than fixed. We assume that sj values are hierarchically 
distributed according to a normal distribution. We take an Empirical Bayesian approach 
in specifying the mean and variance of this distribution, set to the empirical mean and 
empirical variance of the set of maximum titers across sera {max(Hij, . . . , H n j) : j = 
1, . . . , k}. This formulation assumes that particular sera are more reactive in general than 
other sera. 

Additionally, we follow the same logic and assume that some virus isolates are more 
reactive than other virus isolates and include a 'virus effect' Vi representing the general 
level of reactivity across HI assays. With virus reactivity included, observed titers follow 

%~Af(^^-<%, A (10) 

and the vector of virus effects v i for % = 1, . . . , n is estimated in an analogous hierarchical 
fashion, with v normally distributed with mean and variance equal to the empirical mean 
and variance of the set of maximum titers across viruses {max(Hn, . . . , H^) : i = 1, . . . , n}. 



Drift model of antigenic evolution 

As presented, multiple configurations of virus and serum locations X and Y will give the 
same likelihood of an observed data matrix H. An example of this phenomenon is shown 
in Figure [7j In this case, it is impossible to determine from the HI data at hand whether 
the blue and yellow viruses are antigenically similar (Figure Up) or antigenically divergent 
(Figure Up)- This presents an issue of model identifiability, where absolute, as opposed to 
relative, antigenic locations cannot be determined from the observing the serological data 
alone. Thus, in order to achieve more a interpretable model we impose a weak prior on 
global locations. In influenza, it's clear that antigenic distance between strains increases 
with time [6 29 . To capture this, we replace our previous diffuse prior with an informed 
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Figure 7. Schematic antigenic map with three viruses and two sera. (A) Map with 
virus 1 and virus 3 antigenically similar. (B) Map with virus 1 and virus 3 antigenically divergent. 
Virus 1 is shown in blue, virus 2 is shown in red and virus 3 is shown in yellow. Virus isolates are 
represented by filled circles, sera raised against viruses are shown as open circles and map distances 
5ij are shown as solid lines connecting viruses and sera. Sera from virus 1 is compared against 
viruses 1 and 2, while sera from virus 2 is compared against viruses 2 and 3. Configurations (A) 
and (B) represent cartographic models that would give equal likelihoods to a set of serological data 
{Hu,H 2 i,H 2 2, H 32 }. 

prior in which the expected location of viruses and sera increases with date of sampling 
along dimension one, and each virus and serum location follows an independent normal 
distribution centered around this temporal expectation, so that 



x a ~ (3ti+N(0,a 2 x ) 
I/ji ~0tj+-/V(O,aJ), 



(11) 



where t is the difference between the date of the indexed virus or serum and the date 
of the earliest sampled virus or serum, and other dimensions follow Xi m ~ A/"(0, a%.) an d 
Vjm ~ -A/"(0, Cy) for m > 2. Thus, this model assumes that virus and serum locations drift 
in a line across the antigenic map at rate /3. The parameter a x determines the breadth of 
the cloud of virus locations at each point in time, while a y determines the breadth of the 
cloud of serum locations. 



Phylogenetic diffusion model of antigenic evolution 

We simultaneously model antigenic locations and genetic relatedness by assuming that 
virus locations are influenced by evolution following a Brownian motion process 1 10 . To 
do this, we replace the previous prior specifying independent virus locations with a prior 
that incorporates covariance based on shared evolutionary history 



X 



'pt x o N 



K pt n 0, 



+ Evolutionary Brownian Process(cr x ,r) 



(12) 



for P = 2, where a x is the volatility parameter of the Brownian motion and r is a phy- 
logeny specifying tree topology and branch lengths. Thus, viruses which are genetically 
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similar are induced to have prior locations close to one another on the antigenic map. In 
the evolutionary Brownian process, the tips of the phylogeny r correspond to the set of 
virus locations (xi, . . . ,x n ), and the probability of observing tip locations depends on the 
locations of internal nodes (x n+ i, . . . , X2 n _2) and on the location of the root node X2 n _i. 
This process assumes that a virus location x, follows from the location of its parent virus 
xj(j), and with the addition of drift along dimension 1, is distributed as 

x,~(M, 0) / + AA(x /(j) ,d i S) (13) 

for P = 2, where f(i) is a function that maps nodes to parental nodes, d% is the length of the 
branch connecting virus i to parent virus f(i), and XI is a diagonal matrix with diagonal 
elements all equal to a^.. The root virus location X2 n -i is assumed to follow a normal 
distribution with expectation (f3t2 n -i,0)' for P = 2 and variance determined by the 



diffusion volatility a x 10 . The probability of virus locations p(X.\(3,a x ,r) is determined 
through analytical integration across internal states following the methods introduced 
in 10 . This formulation corresponds to a Wiener process with drift, in which the drift 
term j3 only influences the expected states of nodes along the phylogeny, but does not 
influence the covariance structure among these nodes, which remains the same as it does 



in a standard Wiener process [311. This allows the separation in equation 12 between drift 
terms affecting only expectations and the evolutionary Brownian process that includes 
covariance among virus locations xi, . . . ,x n . 

Here, the phylogenetic tree r is estimated using sequence data for viruses 1, . . . , n according 
to well established methods implemented in the software package BEAST [32] . 

Posterior inference 

Top-level priors for 1/V 2 , /3, l/c 2 , and 1/cr 2 are assumed to follow diffuse Gamma(a, b) 
distributions with a = 0.001 and b = 0.001. Under the full model, the posterior probability 
of observing virus and serum locations given immunological data is factored 

p(X, Y|H) oc p(H|X, Y, s, v, <p) p(X|ft a x ,r) p(Y\(3, a y ) p(s, v, <p, (3, a x , a y , r). (14) 

We sample from this posterior distribution using the MCMC procedures implemented 
in the software package BEAST [32] . Metropolis-Hastings proposals include transition 
kernels that translate individual virus and serum locations Xj and y,- and individual virus 
effects Vi and serum effects Sj , and other transition kernels that scale the entire set of virus 
and serum locations X and Y and that scale parameters tp, (3, a x and a y . For the present 
analysis, a two-step approach was taken to sample phylogenies, where a posterior sample 
of phylogenies was gathered using sequence data and then, in the cartographic analysis, 
trees from this set were randomly proposed and accepted following the Metropolis-Hastings 
algorithm |33|. 

Genetic, antigenic and surveillance data 

We compiled an antigenic dataset of hemagglutination inhibition (HI) measurements of 
virus isolates against post-infection ferret sera for influenza A/H3N2 by collecting data 
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from previous publications [mplp)ln4j, NIMR vaccine strain selection reports for 2002 
and 2008-2012 (34H40) and the Feb 2011 VRBPAC report (ill. We queried the Influenza 



Research Database 42 and the EpiFlu Database |43| for HA nucleotide sequences by 



matching strain names, e.g. A/HongKong/1/1968, and only strains for which sequence was 
present was retained. If a strain had multiple sequences in the databases we preferentially 
kept the IRD sequence and preferentially kept the longest sequence in IRD. Sequences 
were aligned using MUSCLE v3.7 under default parameters [44]. This dataset had 2051 
influenza isolates (present as either virus or serum in HI comparisons) dating from 1968 
to 2011. However, the majority of isolates were present from 2002 to 2007. Because 
we are interested in longer-term antigenic evolution, we subsampled the data to have at 
most 20 virus isolates per year, preferentially keeping those isolates with more antigenic 
comparisons. We then kept only those serum isolates that are relatively informative to 
the antigenic placement of viruses, dropping serum isolates that are compared to 4 or 
fewer different virus isolates. This censoring left 402 virus isolates, 519 serum isolates and 
10,059 HI measurements. Each virus isolate was compared to an average of 21.9 serum 
isolates, and each serum isolate was compared to an average of 18.0 virus isolates. 

Antigenic data for influenza A/H1N1 was collected from previous publications |8| |11[[45] - 
58] and NIMR vaccine strain selection reports for 2002-2010 [3lp7l[59H67l. The same 



procedure was followed as was followed for H3N2 to match sequence data and to subsample 
antigenic comparisons. This procedure yielded 115 virus isolates, 77 serum isolates and 
1882 HI measurements over the course of 1977 to 2009. Each virus isolate was compared 
to an average of 10.0 serum isolates, and each serum isolate was compared to an average 
of 16.2 virus isolates. 

Antigenic comparisons for influenza B/Victoria were collated from previous publications 
[2j[8}[68}{75] and vaccine strain selection reports for 2002-2012 [3lH4^[5^[]6^[76}|78) . Here, 



the sequence matching and subsampling procedure yielded 179 virus isolates, 70 serum 
isolates and 2003 HI measurements over the course of 1986 to 2011. Each virus isolate 
was compared to an average of 6.5 serum isolates, and each serum isolate was compared 
to an average of 16.7 virus isolates. 

Antigenic comparisons for influenza B/Yamagata were collected from previous publications 
(2j|8j|68}{75}[79}^5) and vaccine strain selection reports for 2002-2012 [34}|40j[5^}|66j[76}|78 



For B/Yamagata, the matching and subsampling procedure resulted in 174 virus isolates, 
69 serum isolates and 1962 HI measurements over the course of 1987 to 2011. Each virus 
isolate was compared to an average of 6.9 serum isolates, and each serum isolate was 
compared to an average of 17.3 virus isolates. 

Surveillance data was obtained from the Centers of Disease Control and Prevention Flu- 
View Influenza Reports from the yearly summaries of influenza seasons 1997-1998 to 



2010-2011 86 . As an example, one report states "collaborating laboratories in the United 



States tested 195,744 respiratory specimens for influenza viruses, 27,682 (14%) of which 
were positive. Of these, 18,175 (66%) were positive for influenza A viruses, and 9,507 
(34%) were positive for influenza B viruses. Of the 18,175 specimens positive for influenza 
A viruses, 7,631 (42%) were subtyped; 6,762 (87%) of these were seasonal influenza A 
(H1N1) viruses, and 869 (13%) were influenza A (H3N2) viruses." In this case, we esti- 
mate the relative proportion of A/H3N2 of the four clades as 0.66 x 0.13 = 0.09. Similar 
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calculations were performed for A/H1N1, B/Vic and B/Yam. 



Implementation 

Phylogenetic trees were estimated for A/H3N2, A/H1N1, B/Vic and B/Yam using BEAST 
[32] and incorporated the SRD06 nucleotide substitution model [87] , a coalescent demo- 
graphic model with constant effective population size and a strict molecular clock across 
branches. MCMC was run for 60 million steps and trees were sampled every 50,000 steps 
after allowing a burn-in of 10 million steps, yielding a total sample of 2000 trees. These 
trees were treated as a discrete set of possibilities when subsequently sampled in the 
BMDS analysis [33] . However, it would be possible to jointly sample from sequence data 
and serological data using these methods. 

MCMC was used to sample virus locations X, serum locations Y, virus effects v, serum 
effects s, MDS precision 1/f 2 , antigenic drift rate ft, virus location precision l/cr 2 , serum 
location precision 1/cr? and phylogenetic tree r. MCMC chains were run for 500 million 
steps and parameter values sampled every 200,000 steps after a burn-in of 100 million steps, 
yielding a total of 2000 MCMC samples. In all cases when drift parameter ft was included 
the MCMC chain mixed well and arrived at the same estimated posterior distribution 
from different starting points. However, without drift parameter ft, maps for A/H3N2 
showed some degree of metastability, where some chains would converge on one solution 
and other chains would converge on a different solution. We favor models that include 
ft, because its inclusion, in addition to correcting most identifiability issues, yields much 
improved mixing of antigenic locations. 

There is some difficulty summarizing posterior cartographic samples, as sampled virus and 
serum locations represent only relative quantities, and because of this, over the course of 
the MCMC, virus locations may shift. Our prior on virus and serum locations removes 
much of this issue, orienting the antigenic map along dimension 1 and fixing it to begin at 
the origin. However, local isometries are often still a problem. For example, in A/H3N2 the 
HK/68, EN/72 and VI/75 clusters may rotate in relation to other clusters. Consequently, it 
may be difficult to fully align MCMC samples using Procrustes analysis. For the present 
study, we take a simple approach and sample a single MCMC step and visualize the 
antigenic locations at this state (Figures [2j Figure [3]). Then, for specific quantities of 
interest, like rate of antigenic drift and rate of diffusion at different points along the 
phylogeny, we calculate the quantity across MCMC samples to yield an expectation and a 
credible interval. This approach accurately characterizes uncertainty that may be hidden 
in an analysis of a single antigenic map. 



Estimation of diffusion coefficient 

The distance d traveled by a diffusion is not a linear function of time interval t, and 
thus the 'rate' of diffusion cannot be calculated following r = d/t. Assume we have a 
one-dimensional diffusion starting from xq = 0, so that 

X t ~N{0,a 2 t), (15) 
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where a represents the diffusion volatility, and in this case, the expected squared displace- 
ment follows 

m?} = <? 2 t. (i6) 

In moving to a two-dimensional diffusion 

^^^((o)-Ut^))' (17) 

where a x represents volatility in dimension 1, a y represents volatility in dimension 2 and 
p represents correlation, expected squared displacement is equal to the sum of squared 
displacements in each dimension 

E[X 2 + Y 2 ] = a 2 x t + a 2 t. (18) 

Such diffusions are often characterized with diffusion coefficient D. In this case 

D = -^-±, (19) 

so that expected squared displacement follows from D alone 

E[Xf + Y t 2 } = ADt. (20) 



Owing to this relationship, Pybus et al. 20 estimate D from phylogenies in which internal 



nodes have inferred 2D locations via the formula 

in which d% represents the displacement of a branch and t% represents the length of a 
branch. However, uncertainty in di or ti, especially along short branches of the phylogeny, 
may cause df/Ati to have high variance and increase uncertainty in the resulting estimates 
of D. Consequently, we estimate D by looking at the relationship between total squared 
displacement and total time along branches of a phylogeny 

u 

2 



4 = £d; 



i=i 



*E = X}*i 



i=l 



b = A. (22) 



This estimator showed less variance than equation 21 when analyzing tree output from 



the BMDS model and was thus preferred. Additionally, we confirmed through simple 



simulation that estimator 22 shows lower mean squared error than estimator [21] in the 
presence of modest levels of noise on estimates of di and U. When estimating D specific 
to trunk branches or specific to side branches, we calculate d% and is only from a subset 
of phylogeny branches as appropriate. 
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Comparison with previous results 

Here, we attempt to compare antigenic locations inferred by our BMDS model to antigenic 
locations previously inferred by the error minimization methods of Smith et al. |6| , referred 
to here as antigenic cartography by MDS. For this comparison, we use exactly the same HI 
data used to produce the results in p], consisting of 273 virus isolates, 79 serum isolates 
and a total of 4252 HI measurements taken between 1968 and 2003. We begin with a 
BMDS analog of the antigenic model used in [6J, where serum effects are taken as the 
maximum titer of a particular ferret serum and the expected log2 drop in HI titer is 
proportional to Euclidean distance between virus and serum locations. To bring models 
into further alignment, we use a Uniform(— 100, 100) distribution over virus locations and 
serum locations. Unsurprisingly, we find that this BMDS model produces results that are 
strongly congruent with MDS results (Figure [8]). Antigenic cluster locations are consistent 
between methods (Figure p]A-B) and antigenic distances between pairs of viruses are 
consistent between temporally similar and temporally divergent viruses (Figure ^C-E) , 
suggesting that the resulting maps are consistent at both local and global scales. Credible 
intervals of antegenic distances for the BMDS model remain narrow across the temporal 
spectrum (Figure [8p-E) , implying a fair degree of rigidity to the map. 

Smith et al. |6| show that there exist at least two solutions in their assignment of antigenic 
locations, involving the rotation of clusters HK/68, EN/72 and VI/75 (shown in Figure 
S2 of 16)). We observe the same metastable behavior in our analysis; some MCMC chains 
converge on the solution shown in Figure [8^3, while other MCMC chains converge on 
the alternative solution shown in Figure [9j3. The distribution of likelihood values appear 
highly similar between these two solutions, suggesting that they represent global optima. 
The rotation of the HK/68, EN/72 and VI/75 clusters creates a map that bends slightly, 
so that temporally distant viruses appear closer in the rotated solution than in the original 
solution (Figure ^C-E) . In this case, it is clear that the solutions are locally consistent 
between viruses up to ~15 years divergent, even if there some degree of global flexibility. 

As discussed in the main text, the presence of multiple optima with different degrees 
of 2D curvature implies an issue of identifiability; the HI likelihood model alone cannot 
distinguish between these possibilities. Because of this issue, and to more easily estimate 
rates of antigenic drift, we include a model of systematic drift in antigenic location that 
favors linear movement in the antigenic map. We find that including this drift prior on 
antigenic locations removes the problem of identifiability. Antigenic locations produced by 
this model remain locally consistent with MDS results between viruses ~15 years divergent, 
but global comparisons show that this BMDS model has partitioned more variance to the 



first antigenic dimension (Figure 10). We additionally find that including the drift prior 
on antigenic locations often results in greater predictive power, with a slight improvement 
of test error for the A/H1N1, B/Vic and B/Yam datasets (Table [l]). 

Our final BMDS model (model 9, Table [TJ differs from antigenic model used by Smith et 
al. [6] in including temporally- and phylogenetically-informed priors on antigenic locations 
and also in estimating serum and virus effects. Here, we investigate the impact on antigenic 
locations of estimating virus and serum effects in the BMDS model. To isolate this differ- 
ence, we use a Uniform(— 100, 100) prior on antigenic locations. Surprisingly, estimating 
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Figure 8. Comparison of A/H3N2 antigenic locations estimated by Smith et al. J6J 
using MDS and an equivalent BMDS model. (A) MDS antigenic locations, reorientated 
so that the primary dimension lies on the x-axis rather than on the y-axis as in Figure 1 of pi. 
(B) A posterior sample of antigenic locations from an equivalent BMDS model. In (A) and (B), 
viruses are shown as colored circles, with color denoting antigenic cluster inferred by 16], and sera 
are shown as gray points. (C) Antigenic distances between A/Bilthoven/16190/1968 and all other 
viruses determined for both methods. (D) Antigenic distances between A/Fujian/41 1/2002 and all 
other viruses determined for both methods. (E) Antigenic distances between 750 random pairs of 
viruses determined for both methods. In (C), (D) and (E) red points show distances for the MDS 
model and gray bars show the 95% credible interval of distances for the BMDS model, while the 
red dashed line shows a LOESS regression to MDS distances, and the black dashed line shows a 
LOESS regression to the BMDS distances. The BMDS model has a Uniform (-100, 100) prior on 
antigenic locations and serum effects fixed at maximum titer values. 



virus and serum effects results in a more linear antigenic map (Figure 11), resembling the 
appearance of the map incorporating the antigenic drift prior, while preserving local con- 
sistency. We generally observe congruence between MDS and BMDS antigenic locations 
for viruses less than ~10 years divergent (Figure [lip). However, specific viruses may be 
affected, for instance A/Bilthoven/16190/1968 (Figure [lip), which appears more distant 
from all other viruses when serum and virus effects are included. 

In this dataset, viruses 15 or more years divergent always yield threshold titers, and hence, 
their relative locations must be indirectly inferred rather than through direct comparison. 
This may explain why we observe local consistency between models at scales less than 
~15 years, but some degree of global inconsistency. Still, these results suggest that, when 
making local comparisons, such as those used to calculcate year-to-year antigenic drift 
(Figure p]), outcomes are expected to be robust to many model particulars. 
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Figure 9. Comparison of A/H3N2 antigenic locations estimated by Smith et al. [6] 
using MDS and an equivalent BMDS model under an alternative solution. (A) MDS 

antigenic locations, reorientated so that the primary dimension lies on the x-axis rather than on the 
y-axis as in Figure 1 of |6l. (B) A posterior sample of antigenic locations from an equivalent BMDS 
model that has converged on the alternative solution. In (A) and (B), viruses are shown as colored 
circles, with color denoting antigenic cluster inferred by [6], and sera are shown as gray points. (C) 
Antigenic distances between A/Bilthoven/16190/1968 and all other viruses determined for both 
methods. (D) Antigenic distances between A/Fujian/41 1/2002 and all other viruses determined 
for both methods. (E) Antigenic distances between 750 random pairs of viruses determined for 
both methods. In (C), (D) and (E) red points show distances for the MDS model and gray bars 
show the 95% credible interval of distances for the BMDS model, while the red dashed line shows 
a LOESS regression to MDS distances, and the black dashed line shows a LOESS regression to the 
BMDS distances. The BMDS model has a Uniform(— 100, 100) prior on antigenic locations and 
serum effects fixed at maximum titer values. 



Availability 

Source code implementing the cartographic models has been made fully available as part of 
the software package BEAST [32] , and can be downloaded from its Google code repository 
(http://code.google.eom/p/beast-mcmc/). 
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Figure 10. Comparison of A/H3N2 antigenic locations estimated by Smith et al. (6] 
using MDS and an extended BMDS model that includes date-informed priors on 
antigenic locations. (A) MDS antigenic locations, reorientated so that the primary dimension 
lies on the cc-axis rather than on the y-axis as in Figure 1 of [6] . (B) A posterior sample of antigenic 
locations from a BMDS model that includes date-informed priors on antigenic locations. In (A) and 
(B), viruses are shown as colored circles, with color denoting antigenic cluster inferred by |6|, and 
sera are shown as gray points. (C) Antigenic distances between A/Bilthoven/16190/1968 and all 
other viruses determined for both methods. (D) Antigenic distances between A/Fujian/41 1/2002 
and all other viruses determined for both methods. (E) Antigenic distances between 750 random 
pairs of viruses determined for both methods. In (C), (D) and (E) red points show distances for 
the MDS model and gray bars show the 95% credible interval of distances for the BMDS model, 
while the red dashed line shows a LOESS regression to MDS distances, and the black dashed line 
shows a LOESS regression to the BMDS distances. The BMDS model has a date-informed prior 
on antigenic locations and serum effects fixed at maximum titer values. 
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Figure 11. Comparison of A/H3N2 antigenic locations estimated by Smith et al. ro] 
using MDS and an extended BMDS model that estimates serum and virus effects. 

(A) MDS antigenic locations, reorientated so that the primary dimension lies on the x-axis rather 
than on the y-axis as in Figure 1 of (6). (B) A posterior sample of antigenic locations from a 
BMDS model that estimates virus and serum effects. In (A) and (B), viruses are shown as colored 
circles, with color denoting antigenic cluster inferred by Rj], and sera are shown as gray points. (C) 
Antigenic distances between A/Bilthoven/16190/1968 and all other viruses determined for both 
methods. (D) Antigenic distances between A/Fujian/411/2002 and all other viruses determined 
for both methods. (E) Antigenic distances between 750 random pairs of viruses determined for 
both methods. In (C), (D) and (E) red points show distances for the MDS model and gray bars 
show the 95% credible interval of distances for the BMDS model, while the red dashed line shows 
a LOESS regression to MDS distances, and the black dashed line shows a LOESS regression to the 
BMDS distances. The BMDS model has a Uniform(— 100, 100) prior on antigenic locations and 
virus and serum effects estimated in a hierarchical Bayesian fashion. 
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